{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "from scipy.optimize import optimize\n",
    "from scipy.special import expit\n",
    "from matplotlib import pyplot as plt\n",
    "%matplotlib inline\n",
    "from IPython.display import clear_output"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Sample -KL(q || p)\n",
    "def sample_mdkl(log_alpha):\n",
    "    eps = 1 + np.exp(0.5 * log_alpha) * np.random.randn(*log_alpha.shape)\n",
    "    return 0.5 * log_alpha - np.log(np.abs(eps))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Estimate -KL on a grid\n",
    "M = 1000\n",
    "log_alpha_grid = np.linspace(-10, 10, M)\n",
    "mdkl_sampled = np.zeros(M)\n",
    "K = 100000\n",
    "for i in range(100):\n",
    "    print(i)\n",
    "    mdkl_sampled += 0.01 * sample_mdkl(np.tile(log_alpha_grid, K).reshape(K, M)).mean(axis=0)\n",
    "    clear_output()\n",
    "C = sample_mdkl(np.tile([20], 100000000)).mean()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [],
   "source": [
    "def mdkl_approx(log_alpha, k1, k2, k3):\n",
    "    return k1 * expit(k2 + k3 * log_alpha) - 0.5 * np.log1p(np.exp(-log_alpha))\n",
    "\n",
    "def mdkl_approximation_loss(x):\n",
    "    k1, k2, k3 = x\n",
    "    return np.sum((mdkl_sampled - mdkl_approx(log_alpha_grid, k1, k2, k3)) ** 2)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Optimization terminated successfully.\n",
      "         Current function value: 0.009276\n",
      "         Iterations: 188\n",
      "         Function evaluations: 340\n",
      "k1, k2, k3 = 0.636189384063 1.8747510114 1.48840725372\n"
     ]
    }
   ],
   "source": [
    "k1, k2, k3 = optimize.fmin(mdkl_approximation_loss, x0=np.array([0., 0., 0.]))\n",
    "print('k1, k2, k3 =', k1, k2, k3)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "C:\\Users\\dmolc\\Anaconda3\\lib\\site-packages\\matplotlib\\cbook\\deprecation.py:106: MatplotlibDeprecationWarning: Adding an axes using the same arguments as a previous axes currently reuses the earlier instance.  In a future version, a new instance will always be created and returned.  Meanwhile, this warning can be suppressed, and the future behavior ensured, by passing a unique label to each axes instance.\n",
      "  warnings.warn(message, mplDeprecation, stacklevel=1)\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYYAAAEOCAYAAACNY7BQAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMS4wLCBo\ndHRwOi8vbWF0cGxvdGxpYi5vcmcvpW3flQAAIABJREFUeJzt3Xl4VPXZ//H3nY1AWJR9CRqQfZMl\nLAGKWFBww4ogWhcUFZdWq10sarVq6/PUim3tT0XiBm6AjwUXKoqKSKHsEHaQHQLIKkiALGS+vz9m\nEhNMICGZOZPM53VdczFzzplzPnNmmHvOkvuYcw4REZE8UV4HEBGR8KLCICIihagwiIhIISoMIiJS\niAqDiIgUosIgIiKFqDCIiEghKgwiIlKICoOIiBQS43WAs1G3bl2XlJTkdQwRkQpl6dKlB5xz9c40\nXYUsDElJSSxZssTrGCIiFYqZbS/JdNqVJCIihagwiIhIISoMIiJSiAqDiIgUosIgIiKFhEVhMLPB\nZrbBzDaZ2Riv84iIRDLPC4OZRQMvApcB7YAbzKydt6lERCJXOPwdQw9gk3NuC4CZTQauBtZ6mkpE\nfuAc7mQW2VnHyc7KJCc7m5zsLHJysjmZk8XJnGzIPYkvNxuXmwO52biTObjcHJzz4fPl4pzD+XyB\nWy7OBe67XJzP4Xy54Hz+ZblccA7DgfNh5I13+AB8gUsSG/5pTsn6w92C4364b8UML8jnnD+DBR77\nHAaYgRVe4A/zLTivklw22QretfzlOgdRxqkLAqBJv9tofEGHM8+7DMKhMDQBdhZ4nA70PHUiMxsN\njAY477zzQpNMpBJwznHsWAZH9u3k6IFdnPj+ANkZB8k99h3u+HdY1hFisw4Td/IoMb5MYnNPEOvL\noorLJM5lEe+yqEom0eaoAlTx+gVVID5XxDd7Ga1OSomIwlDUmvtRqXXOpQKpAMnJySUoxSKRwTnH\nnn372Lt1Dcf3bMR3aDNx32+n6om91Mg5QG3fIWrZMaoX8/yjVOOo1eB4VHVyouLJjK7JsSrx5EZX\nxRfjv+XGVMPFVIXYeCw6jqiYWKKiY4mKiSMqNo6omDiIisWiYyA6DouOxaJjITqWqKhooqKiiIqK\nxqIMi4omOiqaqOgozKKJjs4bH4VFRWNRUf6bGVgU4B8eZRZ4jvl/TQd+nzsfgZ/XBX/Nm//5UOBf\n//B8VvCr58fDo6MMK7ANEBMVhQNy87cELLD+waKKLgD2o22LglszP57eOf9yowxO+or+mmsfHfwj\nAOFQGNKBpgUeJwK7PcoiEtacL5c9m1eza/1CMnemUf27dTTN3kxjO0LjAtPtow6HY+vxfUIzDlbr\ngS+hIVG1GhFXqxHVzqlPQq061Khdn4Qa51IjOoYanr2iiidUX5ox0SFaUFHL9m7R+RYDLc2sGbAL\nuB74ubeRRMJDZnYOm1fM5cj6r6m2ZyFJx1fSmAwaA9kuhvS4JHbW7cuOuq1IaNSacxNbUzuxFfXj\nq1Pf6/BSYXleGJxzJ83sl8BnQDTwunNujcexRDyzZ98BVv/nA6I3fUqn4wtpb98DkG6N2FCrH7lN\ne1G/TS+SWnWmeZz2+Ev587wwADjnPgE+8TqHiFeysrNZ/vWH5C6fRJdjc7nEssiwBHbU7c2eFoNo\n3GUQiQ3PI9HroBIRwqIwiESqDVu2kf75i3TY8z69OMT3JLCx4RU0TLmeBh1/SrvoWK8jSgRSYRDx\nQNqK5Rz67Bl6H/uC1pbDumrJHOz6JK37DefCuKpex5MIp8IgEkLLV6/h4Cd/5qJjn5Fr0WxpchVN\nBv2atud39DqaSD4VBpEQ2H/oO5a98xj9D0wmyhwbmw6j2dA/0q52E6+jifyICoNIEPl8jjnTJ9Ji\n2Z8ZxH7W1R9MsxHP0K5uktfRRIqlwiASJNt2prPz7Xvpn/U16THnsevK92nb+RKvY4mckQqDSBDM\nnfkvLpj3O3rZYda2uY+2wx7HYuK8jiVSIioMIuXoeFY2c1/7PQP3vsHumEQOj3iXdq16eR1LpFRU\nGETKybf79rL1lZu5NGcha+pfRuvbXyUmvrjWdSLhS4VBpBxsWL+GuMnX0d3t5pvkx2h/5W9O6d4p\nUnGoMIiU0fLFc2k8/SaqWRa7h0ymVbdBXkcSKRMVBpEyWDL7Y1p/dSeZUdXIvvkTzmvexetIImWm\nwiByluZ/+QGd59zJgZgG1LzzY2o1bOZ1JJFyocIgchYWffWhvyjENuTcez6leh39BbNUHioMIqW0\nev5ndJh9J/tjGlLn3s9IqN34zE8SqUCCf/FQkUpk89olnPfpbRyKrkOtu2eoKEilpMIgUkJ7dmwi\n4b0RZFssMSOnUaueLpsjlZMKg0gJHD64nxMThlKdYxwdNpmG57fxOpJI0KgwiJzByZwctqeOIDE3\nnR2XpNKsQ4rXkUSCSoVB5AwWvXo/F2YtZWWnx2jXZ4jXcUSCToVB5DQWfTiO3nvfZVG9a0m+9kGv\n44iEhAqDSDG+WTaHC5c9xtq4TnS9c5zXcURCRoVBpAjfH9xHzY9HccjOodGdU4iJq+J1JJGQUWEQ\nOYXz+djy2khq+w5xZMhrnFtPf6sgkUWFQeQUiyY/Tefj/2Vpq1/TputFXscRCTkVBpECtqTNocuG\nv7O8Wm963fCI13FEPOFpYTCz4Wa2xsx8ZpbsZRaRE98fouqHd3DQzuX8297AovS7SSKT15/81cBQ\nYI7HOSTSOcem12+nru8Aey99idr1GnqdSMQznnZXdc6tAzBdAlE8turTV+l4eBZfJd7Fxb11BTaJ\nbF5vMZSYmY02syVmtmT//v1ex5FK5MCuzSQtfJy10W3pPfLPXscR8VzQC4OZfWFmq4u4XV2a+Tjn\nUp1zyc655Hr16gUrrkQY58tl31u3E+VyqXb9K1SJi/M6kojngr4ryTk3MNjLEDlbae//lS6Zy5nX\n9g/0adnR6zgiYaHC7EoSKW97NqXRds1zLKvSg5Thv/E6jkjY8Pp01WvMLB1IAf5tZp95mUcihy8n\nmxNT7uAE8TS46RWiovUbSSSP12clTQOmeZlBItPKSY/SOWcjc7v+jb5Nk7yOIxJW9DNJIs6u1f+h\nw+ZXmZtwCX2GjPI6jkjYUWGQiHIyMwOm3cU+q02rW1/S39CIFEGFQSLK2jcfpEnuLrb2eZb69ep7\nHUckLKkwSMTYsehjOu1+jy/PGUbvgdd4HUckbKkwSETIPnqIajN+xWYS6XLb37ULSeQ0VBgkImya\neA+1fIfZP/Cf1K5V0+s4ImFNhUEqva1fv027A5/yZf1b6dV3gNdxRMKeCoNUapmHdlF79hjWWgt6\n3/a013FEKgQVBqm8nGPHxDuo4svkxBUvUbNaVa8TiVQIKgxSaX0z40VaHfkvsxJ/Qbfknl7HEakw\nVBikUjq86xuaLPoTy6M78dORf/A6jkiFosIglY7LPcm+t0bhc0bCiFTi42K9jiRSoagwSKWT9n9P\n0ypzFUvbP0yrVm29jiNS4agwSKWSvnY+7dc9z6Kqfel37X1exxGpkFQYpNLIOXEU/nU731ktzh+p\nayyInC39z5FKY+0bv6Dxyd1s+8nfaNCwsddxRCosFQapFDZ+9RYX7vuQ2fVupOcANcgTKQsVBqnw\nDu3eQoOvf8+6qJb0HDXW6zgiFZ4Kg1RouSdPsm/iLUS7XGKGv0aC/rpZpMxUGKRCW/jWH2iTtYpV\nFz5Gy7YXeh1HpFJQYZAKa/W8T+ixbTxLaw6g58/u9TqOSKWhwiAV0v69O2nw+b18G9WQNne8ikXp\noyxSXvS/SSqc3JMn2fPazdRwGZwcNoGEmrW9jiRSqagwSIWzcOIYOmUvZ3Xnx0hqr66pIuVNhUEq\nlKVfTaXXjldZUmsQyT+73+s4IpWSCoNUGNu2bSJp9q9Ij2lKh9GvgpnXkUQqJU8Lg5k9a2brzWyl\nmU0zs3O8zCPh6+jxE3z/1s1Usyzib3yb+ISaXkcSqbS83mL4HOjgnOsEfAM87HEeCUM+n2Pe+Pvp\nlLuWXX3/l/rN9fcKIsHkaWFwzs10zp0MPFwAJHqZR8LTp1NeZPCR91iXeB0tBt7udRyRSs/rLYaC\nRgEzihtpZqPNbImZLdm/f38IY4mX5v93Nhevf5ItVTvS5tYXvI4jEhFigr0AM/sCaFjEqEedcx8G\npnkUOAm8U9x8nHOpQCpAcnKyC0JUCTNbduwk8bM7OR5Vg8aj38NiqngdSSQiBL0wOOcGnm68mY0E\nrgQGOOf0hS8AHD2eycGJN3GhHeLwdR8Sf66uryASKl6flTQY+D0wxDl33MssEj5yfY654++ne24a\nO3v9ifpt+3odSSSieH2M4QWgBvC5maWZ2cse55Ew8PG7L3LZkSlsSBzOBYPVHE8k1IK+K+l0nHMt\nvFy+hJ8vv57FpRufYntCR1rd+qLXcUQiktdbDCL5VnyzhZaz7iIzOkEHm0U85OkWg0ie3QePcPLd\nm2hg35F1w0fEnqODzSJe0RaDeO54Vg4rU++gG2s4NOBv1GzZ2+tIIhFNhUE85fM5/v3K4wzOmsm2\ndvfQ6Ce3eB1JJOKpMIinPnp/AkP3j2NrvYtJGvY/XscREVQYxEOz/jOHn655mD1VLyDpjrdBl+cU\nCQv6nyieSNuwmRZf3I4vOp56d07FqlT3OpKIBOisJAm5nfsPkzvJfwZS5g0fUaXO+V5HEpECtMUg\nIfX9iWxWpd5JN9Zy+JK/U0tnIImEHRUGCZmcXB8fjX+My3Nmkt7hHhr0udnrSCJSBBUGCQnnHO++\n/Qo3fDee9IYDSByqM5BEwpUKg4TEB598wrAtj7G/eisSR72lM5BEwpj+d0rQzVmSRu9F95IVW4v6\noz+EuASvI4nIaeisJAmqtVt30vDjm6gelU30rdOJqtXI60gicgbaYpCg+fbQUY6+eSPNbTfZ104g\nPrGj15FEpATKVBjMbF55BZHK5VhmDivGj6KnW8G+i/7CuR0GeR1JREqorFsM6o0sP5Lrc3w2/iEG\nZc1kW7t7aXzxaK8jiUgpnPEYg5m9AKwEVgGrnXNHC4x2wQomFddHb/2Dod+9zpaGl9N8uE5LFalo\nSnLweSXQCbgR6GBm3+MvEqvwX69ZJN/MGVO5fMuf2VajM83vmABmXkcSkVIqSWFIcc7dlvfAzBLx\nF4qOwKfBCiYVz5KlC+m+4D4OxTYi8e6poEtzilRIJTnGEG1mT+U9cM6lO+c+AcYB8UFLJhXKlm3b\naPDRzVhUNDXv+ICY6nW8jiQiZ6kkhWEUkGJmd+QNMLMLgUXAumAFk4rjwHeHOT5xOPXsO7Kve5eE\nhi28jiQiZXDGXUnOuZNmNhSYbWa7gSbAo8Adzrkvgh1Qwltmdg7fvPxzevk2smPASyS17et1JBEp\no5KclfQS/gPQfwXewH/Quadzbm+Qs0mY8/kcc8f9goFZ81h/4e9p0+/nXkcSkXJQkoPPafxwsDkW\naAW8YmargFXOuclBzCdhbNbb/8PA76awqvF1dLzmYa/jiEg5KcmupNSCj085K+kKQIUhAs375F0u\n3vwsa2v2psPtL+m0VJFK5Gya6HVxzn0MfFLWhZvZn4CrAR+wD7jVObe7rPOV4Fq1eA6dFz7A9rgL\naHHPFCw61utIIlKOzqYlxtPluPxnnXOdnHOdgenA4+U4bwmCnVu/ocG/byEjqgZ17/yAuGo1vY4k\nIuXsbApDue0zcM59X+BhAmqxEdYOf3eAnLeGUZVMTo6YQs36Tb2OJCJBcDa7ksr1y9vMngZuAY4A\nF59mutHAaIDzzjuvPCNICWRnZbHj5eG0zU1n86UTaNMm2etIIhIkQb8eg5l9YWari7hdDeCce9Q5\n1xR4B/hlcfNxzqU655Kdc8n16tULdmwpwPl8LB93G52ylrGqyxO06TPE60giEkRBv4Kbc25gCSd9\nF/g38McgxpGzMH/io/Q+/G8WNh1Fz5/d73UcEQmys9liKLc/bDOzlgUeDgHWl9e8pXws/ng8vbe/\nxNKaA+lx23NexxGRECj1FoNz7pJyXP5fzKw1/tNVtwN3l+O8pYxW/3cGFy55hHVVOtLx3rexKF0J\nViQSBH1X0uk45671cvlSvO0b0kiceQffRjegyV1TiYuv6nUkEQkR/QSUHzmwdxcxk6/DRxQxN/+L\nmnXqex1JREJIhUEKOXEsgwOvDKWO7xAHr5pI42ZtvY4kIiGmwiD5cnNzWfvS9bTK2cD63s/RsttP\nvY4kIh5QYZB8C1Pvo9ux/7C41a/pPGik13FExCMqDALA/MnP0HvvOyyqdy09f/6Y13FExEMqDMLy\nL6fQY93/sqJqL7rdNV4ttEUinApDhNuYNo9Wc+5jW0xzWv5iCtExaqEtEulUGCLYnh0bOeeDGzlq\nNah1x1SqVT/H60giEgZUGCLUkcOHODHhWuLJJOu6SdRtlOR1JBEJEyoMESg7K4vt44bRNDedHQPG\nc367Hl5HEpEwosIQYZzPx7Jxo+iUtZQVnZ+g/U+u9jqSiIQZFYYIM//NP9Dr8HQWJI4i+Rq10BaR\nH1NhiCCLPx5P720vsrTmQHqOUgttESmaCkOEWDPf30J7bZxaaIvI6enbIQJs35BGk8/8LbQT1UJb\nRM5AhaGSO7hPLbRFpHRUGCqxE8cy2J96LXV8hzigFtoiUkIqDJWUv4X2DbTKWc/63s/RSi20RaSE\nVBgqKX8L7TlqoS0ipabCUAktmPJXfwvtukPVQltESk2FoZJJ+3IK3df+Dyuq9qTb3alqoS0ipabC\nUIlsTJtHy/wW2u+phbaInBUVhkri252bqPXBTWqhLSJlpsJQCXx/5BDH3xhKVU6QqRbaIlJGKgwV\nXHZWFtte8rfQ3j5gPElqoS0iZaTCUIE5n4/lgRbaaRf+kQ5qoS0i5SAsCoOZ/dbMnJnV9TpLRbLg\nzT/QM9BCu/vQX3kdR0QqCc8Lg5k1BS4BdnidpSJZ/HEqKWqhLSJB4HlhAP4OPAQ4r4NUFP4W2g+z\nLq6DWmiLSLnz9BvFzIYAu5xzK7zMUZFs/2ZFfgvtJndNUwttESl3McFegJl9ATQsYtSjwCPApSWc\nz2hgNMB5551XbvkqkoP7dhE9SS20RSS4gl4YnHMDixpuZh2BZsAK87dtSASWmVkP59y3RcwnFUgF\nSE5OjrjdTpnH/S20k3wH2THkPVqphbaIBEnQC0NxnHOrgPyfvGa2DUh2zh3wKlO48uXmsubFG+iS\ns54VKc/TRS20RSSIdNSyAlj4Sl4L7QfpMlgttEUkuDzbYjiVcy7J6wzhaMF7fyXlW38L7R43qIW2\niASfthjCWNqXU+i+5ocW2jotVURCQd80YWrTinm0UgttEfGACkMY+nbnJmpOu4nvrQa1blcLbREJ\nLRWGMONvoX0tVV2ghXbjJK8jiUiEUWEIIznZeS20d7J9oFpoi4g3VBjChPP5WKYW2iISBlQYwsSC\ntx6j53fTWZB4m1poi4inVBjCwJLpr5Cy9YVAC+2/eR1HRCKcCoPH1sz/lE6Lx7AurgMd7n1Lf6sg\nIp7Tt5CHdmxcQZPPbs9voV0lvprXkUREVBi8cmjfLqLe9bfQjr5JLbRFJHyoMHgg83gG+1Kvpa7v\nIPuvnECT5mqhLSLhQ4UhxPJaaLfKWc+6lLG0Th7gdSQRkUJUGEJs4av3+1tot3yALoNv9TqOiMiP\nqDCE0IL3/krKnrf9LbR//rjXcUREiqTCECIrZqmFtohUDPp2CoFNK+bR8uv72BrTnBb3qoW2iIQ3\nFYYgK9hC+5zbp5JQQy20RSS8qTAEUcEW2ieGq4W2iFQMKgxBkpOdxbZxw/0ttAe8TLP2aqEtIhWD\nCkMQOJ+P5eNG0SlzCWkXPk6Hfj/zOpKISImpMATBgrceo0d+C+0HvI4jIlIqKgzlLK+F9hK10BaR\nCkqFoRytXfAZnRaPYW1cBzqqhbaIVFAxXgeoLHZsXEHjT0cFWmhPVQttqZRycnJIT08nMzPT6yhy\nGvHx8SQmJhIbe3Z/M6XCUA5ObaFdq04DryOJBEV6ejo1atQgKSkJM/M6jhTBOcfBgwdJT0+nWbNm\nZzUPT/d1mNkTZrbLzNICt8u9zHM21EJbIklmZiZ16tRRUQhjZkadOnXKtFUXDlsMf3fOjfU6xNnw\n5eay5qWf0yVnPWkp/6CrWmhLBFBRCH9lfY90dLQMFr56P90yvmZxywfoqhbaIiHz9NNP0759ezp1\n6kTnzp1ZuHBh0JbVv39/lixZUuLpZ8+ezZVXXhm0PKEQDlsMvzSzW4AlwG+cc995HagkFrz3LCl7\n3mahWmiLhNT8+fOZPn06y5Yto0qVKhw4cIDs7GyvY1UqQd9iMLMvzGx1EbergXHABUBnYA/w3Gnm\nM9rMlpjZkv379wc79mmtmPUe3dc87W+hfdd4nZYqEkJ79uyhbt26VKlSBYC6devSuHFjnnrqKbp3\n706HDh0YPXo0zjnA/4v/wQcfpF+/frRt25bFixczdOhQWrZsyR/+8AcAtm3bRps2bRg5ciSdOnVi\n2LBhHD9+/EfLnjlzJikpKXTt2pXhw4eTkZEBwKeffkqbNm3o27cvU6dODdGaCB7LW3leM7MkYLpz\nrsOZpk1OTnal2bQrT5tW/pdG/7qGPTFNaPTALHVLlYiybt062rb1n2Dx5MdrWLv7+3Kdf7vGNfnj\nVe1PO01GRgZ9+/bl+PHjDBw4kBEjRnDRRRdx6NAhateuDcDNN9/Mddddx1VXXUX//v3p2bMnzzzz\nDM8//zzPPPMMS5cupXbt2lxwwQWsWLGCo0eP0qxZM+bOnUufPn0YNWoU7dq147e//S39+/dn7Nix\nJCUlMXToUGbMmEFCQgLPPPMMWVlZPPTQQ7Rs2ZJZs2bRokULRowYwfHjx5k+fXq5rpvSKvhe5TGz\npc655DM91+uzkhoVeHgNsNqrLCXx7c7N1Jx6I0etulpoi3ikevXqLF26lNTUVOrVq8eIESOYMGEC\nX331FT179qRjx47MmjWLNWvW5D9nyJAhAHTs2JH27dvTqFEjqlSpQvPmzdm5cycATZs2pU+fPgDc\ndNNNzJ07t9ByFyxYwNq1a+nTpw+dO3dm4sSJbN++nfXr19OsWTNatmyJmXHTTTeFaE0Ej9fHGP5q\nZp0BB2wD7vI2TvGOHjnE8TeGUt+dYP91H9JMLbQlwp3pl30wRUdH079/f/r370/Hjh0ZP348K1eu\nZMmSJTRt2pQnnnii0OmaebudoqKi8u/nPT558iTw4zN5Tn3snOOSSy5h0qRJhYanpaVVujO1PN1i\ncM7d7Jzr6Jzr5Jwb4pzb42We4uRkZ7E10EJ724CXada+p9eRRCLWhg0b2LhxY/7jtLQ0WrduDfiP\nN2RkZPD++++Xer47duxg/vz5AEyaNIm+ffsWGt+rVy/mzZvHpk2bADh+/DjffPMNbdq0YevWrWze\nvDn/uRWd11sMYS+vhXaPzCUs7vQk3dVCW8RTGRkZ3HfffRw+fJiYmBhatGhBamoq55xzDh07diQp\nKYnu3buXer5t27Zl4sSJ3HXXXbRs2ZJ77rmn0Ph69eoxYcIEbrjhBrKysgD485//TKtWrUhNTeWK\nK66gbt269O3bl9Wrw3qv+BmFzcHn0gjlwef5Ex8lZesLLEi8jV53/CMkyxQJV0Ud0KwMtm3bxpVX\nXlnhv9ALqrAHn8Pdkn+rhbaIRB4VhmKsW/gZHRc9zNpYtdAWqeySkpIq1dZCWenbrgg7Nq6k0YxR\n7IuqR5O71UJbRCKLCsMp/C20h+OIIvpmtdAWkcijwlBAwRba+658g8bN23kdSUQk5FQYAny5uax9\n6ee0ylnP2pSxtE4e6HUkERFPqDAELHz1frpmfM0itdAWCXvTpk3DzFi/fr1nGR5//HG++OKLMs/n\n8OHDvPTSS6V+3hNPPMHYscG5lI0KA7DwvbH5LbR7qoW2SNjL+8vkyZMnl8v88tpilMZTTz3FwIFl\n37NwtoUhmCK+MKz46v/ophbaIhVGRkYG8+bN47XXXssvDLNnz6Zfv35cc801tGvXjrvvvhufzwf4\nm+795je/oWvXrgwYMIC8tv39+/fnkUce4aKLLuL5559n+/btDBgwgE6dOjFgwAB27NgBwNVXX82b\nb74JwPjx47nxxhsBuPXWW/NbbyQlJfHII4+QkpJCcnIyy5YtY9CgQVxwwQW8/PLL+bkHDBhA165d\n6dixIx9++CEAY8aMYfPmzXTu3Jnf/e53ADz77LN0796dTp068cc//jH/tT/99NO0bt2agQMHsmHD\nhqCt44huibF55X9pMfuXbItpRot73yMmNs7rSCIVx4wx8O2q8p1nw45w2V9OO8kHH3zA4MGDadWq\nFbVr12bZsmUALFq0iLVr13L++eczePBgpk6dyrBhwzh27Bhdu3blueee46mnnuLJJ5/khRdeAPy/\n1r/++msArrrqKm655RZGjhzJ66+/zv33388HH3xAamoqffr0oVmzZjz33HMsWLCgyFxNmzZl/vz5\nPPjgg9x6663MmzePzMxM2rdvz9133018fDzTpk2jZs2aHDhwgF69ejFkyBD+8pe/sHr1atLS0gD/\nNR82btzIokWLcM4xZMgQ5syZQ0JCApMnT2b58uWcPHmSrl270q1bt/Ja84VEbGHYm76ZGlNvJMMS\nqKUW2iIVxqRJk3jggQcAuP7665k0aRJXXHEFPXr0oHnz5gDccMMNzJ07l2HDhhEVFcWIESMAfzvt\noUOH5s8rbzj4rwyXd5Gdm2++mYceegiABg0a8NRTT3HxxRczbdq0/Gs+nKpga++MjAxq1KhBjRo1\niI+P5/DhwyQkJPDII48wZ84coqKi2LVrF3v37v3RfGbOnMnMmTPp0qUL4N/S2LhxI0ePHuWaa66h\nWrVqhZYXDBFZGI4eOcSx19VCW6RMzvDLPhgOHjzIrFmzWL16NWZGbm4uZsbll19+xrbZRQ1PSEgo\ndlkFp1u1ahV16tRh9+7dxU5/ptbe77zzDvv372fp0qXExsaSlJRUqDV4HuccDz/8MHfdVfgqBP/4\nxz9C1t474naoF2yhvXXAOLXQFqlA3n//fW655Ra2b9/Otm3b2LlzZ/6V1xYtWsTWrVvx+XxMmTIl\nv222z+fLPxbw7rvv/qiddp4LI+eAAAAJjUlEQVTevXvnH7N455138qdbtGgRM2bMYPny5YwdO5at\nW7eeVfYjR45Qv359YmNj+eqrr9i+fTsANWrU4OjRo/nTDRo0iNdffz3/sqG7du1i37599OvXj2nT\npnHixAmOHj3Kxx9/fFY5SiKithicz8fyl2+nR+YSFnV6kh79rvE6koiUwqRJkxgzZkyhYddeey3j\nxo0jJSWFMWPGsGrVqvwD0eDfKlizZg3dunWjVq1aTJkypch5//Of/2TUqFE8++yz1KtXjzfeeIOs\nrCzuvPNO3njjDRo3bsxzzz3HqFGjmDVrVqmz33jjjVx11VUkJyfTuXNn2rRpA0CdOnXo06cPHTp0\n4LLLLuPZZ59l3bp1pKSkAP6D52+//TZdu3ZlxIgRdO7cmfPPP5+f/OQnpc5QUhHVdnv+m4+RsuWf\nzG9yGyl3qoW2SGmFa9vt2bNnM3bs2CKvs1y9evX8X9+RRG23SyiuzvksOudyet2uFtoiIsWJqF1J\n3a64A7jD6xgiUs7yrv9clEjcWiiriNpiEBGRM1NhEJFSqYjHJSNNWd8jFQYRKbH4+HgOHjyo4hDG\nnHMcPHiQ+Pj4s55HRB1jEJGySUxMJD09Pb/fkISn+Ph4EhMTz/r5KgwiUmKxsbE0a9bM6xgSZNqV\nJCIihagwiIhIISoMIiJSSIVsiWFm+4HtZ/n0usCBcoxTXpSrdJSrdJSrdMI1F5Qt2/nOuXpnmqhC\nFoayMLMlJekVEmrKVTrKVTrKVTrhmgtCk027kkREpBAVBhERKSQSC0Oq1wGKoVylo1ylo1ylE665\nIATZIu4Yg4iInF4kbjGIiMhpVMrCYGbDzWyNmfnMLPmUcQ+b2SYz22Bmg4p5fjMzW2hmG81sipnF\nBSHjFDNLC9y2mVlaMdNtM7NVgelKf9m60ud6wsx2Fch2eTHTDQ6sw01mNqaoaco517Nmtt7MVprZ\nNDM7p5jpQrK+zvT6zaxK4D3eFPgsJQUrS4FlNjWzr8xsXeDz/6sipulvZkcKvL+PBztXYLmnfV/M\n75+B9bXSzLqGIFPrAushzcy+N7MHTpkmZOvLzF43s31mtrrAsNpm9nngu+hzMzu3mOeODEyz0cxG\nljmMc67S3YC2QGtgNpBcYHg7YAVQBWgGbAaii3j+e8D1gfsvA/cEOe9zwOPFjNsG1A3hunsC+O0Z\npokOrLvmQFxgnbYLcq5LgZjA/WeAZ7xaXyV5/cC9wMuB+9cDU0Lw3jUCugbu1wC+KSJXf2B6qD5P\nJX1fgMuBGYABvYCFIc4XDXyL/zx/T9YX0A/oCqwuMOyvwJjA/TFFfe6B2sCWwL/nBu6fW5YslXKL\nwTm3zjm3oYhRVwOTnXNZzrmtwCagR8EJzMyAnwLvBwZNBH4WrKyB5V0HTArWMoKgB7DJObfFOZcN\nTMa/boPGOTfTOXcy8HABcPatI8uuJK//avyfHfB/lgYE3uugcc7tcc4tC9w/CqwDmgRzmeXoauBN\n57cAOMfMGoVw+QOAzc65s/3D2TJzzs0BDp0yuODnqLjvokHA5865Q86574DPgcFlyVIpC8NpNAF2\nFniczo//49QBDhf4EipqmvL0E2Cvc25jMeMdMNPMlprZ6CDmKOiXgc3514vZdC3JegymUfh/XRYl\nFOurJK8/f5rAZ+kI/s9WSAR2XXUBFhYxOsXMVpjZDDNrH6JIZ3pfvP5MXU/xP868WF95Gjjn9oC/\n8AP1i5im3NddhW27bWZfAA2LGPWoc+7D4p5WxLBTT8sqyTQlUsKMN3D6rYU+zrndZlYf+NzM1gd+\nWZy10+UCxgF/wv+a/4R/N9eoU2dRxHPLfHpbSdaXmT0KnATeKWY25b6+iopaxLCgfY5Ky8yqA/8C\nHnDOfX/K6GX4d5dkBI4ffQC0DEGsM70vXq6vOGAI8HARo71aX6VR7uuuwhYG59zAs3haOtC0wONE\nYPcp0xzAvxkbE/ilV9Q05ZLRzGKAoUC308xjd+DffWY2Df9ujDJ90ZV03ZnZK8D0IkaVZD2We67A\nQbUrgQEusHO1iHmU+/oqQklef9406YH3uRY/3k1Q7swsFn9ReMc5N/XU8QULhXPuEzN7yczqOueC\n2heoBO9LUD5TJXQZsMw5t/fUEV6trwL2mlkj59yewK61fUVMk47/WEieRPzHV89apO1K+gi4PnDG\nSDP8lX9RwQkCXzhfAcMCg0YCxW2BlNVAYL1zLr2okWaWYGY18u7jPwC7uqhpy8sp+3WvKWZ5i4GW\n5j97Kw7/ZvhHQc41GPg9MMQ5d7yYaUK1vkry+j/C/9kB/2dpVnHFrLwEjmG8Bqxzzv2tmGka5h3r\nMLMe+L8DDgY5V0nel4+AWwJnJ/UCjuTtQgmBYrfavVhfpyj4OSruu+gz4FIzOzew6/fSwLCzF4qj\n7aG+4f9CSweygL3AZwXGPYr/jJINwGUFhn8CNA7cb46/YGwC/g+oEqScE4C7TxnWGPikQI4Vgdsa\n/LtUgr3u3gJWASsDH8pGp+YKPL4c/1kvm0OUaxP+/ahpgdvLp+YK5foq6vUDT+EvXADxgc/OpsBn\nqXkI1lFf/LsQVhZYT5cDd+d9zoBfBtbNCvwH8XuHIFeR78spuQx4MbA+V1HgbMIgZ6uG/4u+VoFh\nnqwv/MVpD5AT+P66Hf9xqS+BjYF/awemTQZeLfDcUYHP2ibgtrJm0V8+i4hIIZG2K0lERM5AhUFE\nRApRYRARkUJUGEREpBAVBhERKUSFQUREClFhEBGRQlQYRErIzDK8ziASCioMIiJSiAqDSCmZ2a/N\nbHXg9kCB4Y+Z/ypzn5vZJDP7bYFx7c3sCzP7JjDd/zOz7t68ApHTq7DdVUW8YGbdgNuAnvj7+yw0\ns6/xXwHsWvzXQIjB3655aeA5eX2ThuO/utZ6YKlzbnHIX4BICagwiJROX2Cac+4YgJlNxX+xpSjg\nQ+fcicDwjws8ZyCw3Dm3JjAuDv91LkTCknYliZROcZfnPN1lO7vg34LAzBoDGc65eeUdTKS8qDCI\nlM4c4GdmVi1wbYFrgP8Ac4GrzCw+cAW1Kwo8J4sfrlH9v0BcKAOLlJZ2JYmUgnNumZlN4IcLPL3q\nnFsOYGYf4e/bvx1Ygv86zwDvAh+a2QZgPFDFzP7hnHsAkTCk6zGIlBMzq+781wauhn/LYrRzbpnX\nuURKS1sMIuUn1cza4b9620QVBamotMUgIiKF6OCziIgUosIgIiKFqDCIiEghKgwiIlKICoOIiBSi\nwiAiIoWoMIiISCEqDCIiUsj/B+Yqxc2Se2d0AAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x28a3a1a2550>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.plot(log_alpha_grid, mdkl_sampled - C, label='Sampled')\n",
    "plt.plot(log_alpha_grid, mdkl_approx(log_alpha_grid, k1, k2, k3) - C, label='Approximated')\n",
    "_ = plt.axes().set_xlabel(r'$\\log\\alpha$') \n",
    "_ = plt.axes().set_ylabel(r'$-KL$') \n",
    "_ = plt.legend(loc=4)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "C:\\Users\\dmolc\\Anaconda3\\lib\\site-packages\\matplotlib\\cbook\\deprecation.py:106: MatplotlibDeprecationWarning: Adding an axes using the same arguments as a previous axes currently reuses the earlier instance.  In a future version, a new instance will always be created and returned.  Meanwhile, this warning can be suppressed, and the future behavior ensured, by passing a unique label to each axes instance.\n",
      "  warnings.warn(message, mplDeprecation, stacklevel=1)\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAZoAAAEOCAYAAACw8dE2AAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMS4wLCBo\ndHRwOi8vbWF0cGxvdGxpYi5vcmcvpW3flQAAIABJREFUeJzt3Xd8m+W1wPHf8V6xHTvOspM4hAwS\nViYbwt6EMkNLm9JB6YXLKu2FS8ssvVAKXUALlFVKmR2EEkaAsAJkkgSSEOIMiDOdZTtxvM/9433l\nyLJkvba17Jzv56OPpVfveCzLOnrWeURVMcYYY6IlKd4FMMYY07NZoDHGGBNVFmiMMcZElQUaY4wx\nUWWBxhhjTFRZoDHGGBNVFmiMMcZElQUaY4wxUWWBxhhjTFSlxLsAiaBPnz5aWloa72IYY0y3smDB\ngq2qWhRuPws0QGlpKfPnz493MYwxplsRka+87GdNZ8YYY6LKAo0xxpioskBjjDEmqizQGGOMiSoL\nNMYYY6LKAo0xxpioskBjjDEmqizQmB7ji01VfLCyIt7FMMYEsAmbpttraGrmnS+28KOnFwDw+e2n\nkpNub21jEoX9N5puq7lZqWlo4sBb32i1/fXPN3HB+JI4lcoYE8iazky3dc3zi1oFmfysVACWbqiM\nV5GMMUFYoDHd1iuLN7R6/OkvTmZYUTZPzF5L2ZbqOJXKGBPIAo3pliqq61o9TktJQkSorm0E4OTf\nvk9zs8ajaMaYANZHY7qd2oYmzvrjBwA8+8PDGdEvh+QkAeDmMw/gmucWoQqzV23lmOFhM5gbY6LM\najSm25m/dgebq5wazeiBuRTmpJOflQbAlEOL+eSmEwH49mNzqaxpiFs5jTEOCzSm2/nnp+UA/Oi4\n/cjLTG3zfL/c9Jb7H6/eFrNyGWOCs0BjupXKPQ38Z8lGvnXYYG46/YCg+4gI549zhjc/+dGaWBbP\nGBOEBRrTrby8aD31jc1cPHFQu/v99NSRAHyyenssimWMaUdcA42InCYiK0SkTERuDPJ8uog87z4/\nR0RK/Z67yd2+QkRO9dt+nYgsFZHPReRZEcmIzW9jYuH1zzcxqn8vDi7Jb3c/35waY0z8xS3QiEgy\n8CBwOjAauERERgfs9n1gh6ruD/wWuMc9djQwFRgDnAY8JCLJIlIMXA1MUNUDgWR3P9MDVNU2sKS8\nknFDeofdNyM1mfQU5+1tw5yNia941mgmAWWqulpV64HngCkB+0wBnnLvvwScKCLibn9OVetUdQ1Q\n5p4PnCHbmSKSAmQBGzA9wl8/Wsuuukamhmk28/E1n23bXR/NYhljwohnoCkG1vk9Lne3Bd1HVRuB\nSqAw1LGquh74DfA1sBGoVNU3o1J6E3PLN1UztE922GYzn5H9ewHw9vLN0SyWMSaMeAYaCbItsI0j\n1D5Bt4tIb5zazlBgIJAtIpcGvbjI5SIyX0TmV1RYavlE9+HKrby6ZCNpyd7fskcO60OSQPmOPVEs\nmTEmnHgGmnLAvw2khLbNXC37uE1hecD2do49CVijqhWq2gD8Ezgy2MVV9RFVnaCqE4qKbPZ4ovvX\np+sBGFOc6/mY5CShMCedrbvqwu9sjImaeAaaecBwERkqImk4nfbTA/aZDkxz718AvKOq6m6f6o5K\nGwoMB+biNJkdLiJZbl/OicDyGPwuJsq2VNeSm5HCrWeP6dBxRTnpbK6qjVKpjDFexC3Xmao2ishV\nwBs4o8MeV9WlInIHMF9VpwOPAU+LSBlOTWaqe+xSEXkBWAY0AleqahMwR0ReAha62z8FHon172Yi\n6+3lm/lg5Va+f/TQoJkA2jOkMIvlG6uiVDJjjBdxTaqpqjOAGQHbbvG7XwtcGOLYu4C7gmy/Fbg1\nsiU18fTJ6m2kJAk3nj6qw8eO7N+L1z7fxPy125lQWhCF0hljwrHMACah1TU28egHa2hsVlI7MBDA\n55xDBgLwxEdrI1wyY4xXFmhMQlu+0VnA7OIJ3ubOBNqvKIeJpb2pqLIBAcbEiwUak9D+/el60pKT\n+MmpIzp9jn65GVTYyDNj4sYCjUlo731Zwfghvenbq/Mp64p7Z/LVtt1s2GnzaYyJBws0JmEt/HoH\na7buZuxgb5kAQjlhZF+aFVZsqo5QyYwxHWGBxiSsxet2AvCdI0q7dJ4+vZyF0KpqbbVNY+LBAo1J\nWCs2VVOQndZqxczOyM1w5t5U7rFAY0w8WKAxCamxqZl3V1Qwsl8vnCQPneeb5FllgcaYuLBAYxLS\nvW+uYFNVbUuzV1ekpSSRnZZMRbWNPDMmHizQmIT0n8UbARjlpvrvqgOL81j49c6InMsY0zEWaEzC\naWpW6hqbKc7P5EfH7heRcx4wIJe1W3dH5FzGmI6xQGMSztw129m6q46bzhhFSifSzgRTnJ9JdV2j\nDQgwJg4s0JiE89rnG8lMTeaEUX0jds6B+ZkANmnTmDiwQGMSzry1O5hQ2pustMglFx+Y72QWsEBj\nTOxZoDEJpbahiS83V3NwSV5Ez1vc26nRfL7e1qYxJtYs0JiE8vGqbTQ1KwcVdy3tTKCinHT265PN\nrBVbInpeY0x4FmhMQvn3ovXkZqRw9PA+ET2viDB6YK5N2jQmDizQmISyYlM144b0Jic98ou/9spI\npaq2MeLnNca0zwKNSRj1jc18samaEf0iM0kzUG5miiXWNCYO4hpoROQ0EVkhImUicmOQ59NF5Hn3\n+TkiUur33E3u9hUicqrf9nwReUlEvhCR5SJyRGx+G9NV9725AiB6gSYjlfrGZuoam6JyfmNMcHEL\nNCKSDDwInA6MBi4RkdEBu30f2KGq+wO/Be5xjx0NTAXGAKcBD7nnA/g98LqqjgIOAZZH+3cxkfHJ\n6m0AnHXwgKicv3dWGgDbdtVH5fzGmODCBhoROUpEZorIlyKyWkTWiMjqCFx7ElCmqqtVtR54DpgS\nsM8U4Cn3/kvAieKk8p0CPKeqdaq6BigDJolILnAs8BiAqtarqiW46gZq6htZtrGKK44bRkZqcvgD\nOmFIYRYAa7dZKhpjYslLj+tjwHXAAiCSbQ7FwDq/x+XAYaH2UdVGEakECt3tnwQcWwzsASqAJ0Tk\nELfM16iqfbIkuPlrd9DQpBwxrDBq1xhc4ASaddtrYFjULmOMCeCl6axSVV9T1S2qus13i8C1gy0y\noh73CbU9BRgH/ElVxwK7gTZ9PwAicrmIzBeR+RUVFd5LbaLio1XbSEkSJpb2jto1itwlB2y5AGNi\ny0ugmSUi94rIESIyzneLwLXLgUF+j0uADaH2EZEUIA/Y3s6x5UC5qs5xt7+EE3jaUNVHVHWCqk4o\nKirq4q9iumrh1zs4sDgvomlnAmWkJpObkWKBxpgY8/Jf7WvOmuC3TYETunjtecBwERkKrMfp3P9m\nwD7TgWnAx8AFwDuqqiIyHfi7iNwPDASGA3NVtUlE1onISFVdAZwILOtiOU2ULd1Qydw12zkzSoMA\n/BX1SqdilwUaY2IpbKBR1eOjcWG3z+Uq4A0gGXhcVZeKyB3AfFWdjtM/9LSIlOHUZKa6xy4VkRdw\ngkgjcKWq+vqP/ht4RkTSgNXAZdEov4mc5+c5XXVnHRT9QNMnJ52t1TbqzJhYChtoRCQPuBVnNBfA\ne8AdqlrZ1Yur6gxgRsC2W/zu1wIXhjj2LuCuINsX0br2ZRLcvLU7OHr/Ppweg0BT1CudpRsssaYx\nseSlj+ZxoBq4yL1VAU9Es1Bm31G5p4EvNlUxsbQgJtcr6pXOlqpaVAPHnRhjosVLoBmmqre6811W\nq+rtQGTW1zX7vIVf70CVqI4283dwSR6765tY8NWOmFzPGOMt0OwRkaN9D0TkKJz5KsZ02ZtLN5GS\nJBw6OLLLAoRy4EBnnZuNlbUxuZ4xxtuosx8DT7l9NYLTKf/daBbK7Bv21Dfxr0/X842xxVEd1uyv\nV0YqANWWxdmYmPEy6mwRcIib3gVVtZ5UExGzy7ZS29DMOYcOjNk1e2U4b/lqy+JsTMyEDDQicqmq\n/k1Erg/YDoCq3h/lspkebuayzfRKT+GwodFLOxMoKy2Z5CSxGo0xMdRejSbb/RksZ7sN2TFdoqq8\nv7KCo4f3IS0ldknERYSc9BSr0RgTQyEDjao+7N59S1Vn+z/nDggwptPWbqthY2UtVx4f2SWbveiV\nkWI1GmNiyMtXyT963GaMZx+t2grAkVHM1hxKT1vS+cFZZVz2xFxKb3yVFZuq410cY9por4/mCOBI\noCignyYXJ2WMMZ02c9lm+udmMLRPdvidI6xXegobdvacEfr3vrGi5f4HKysY2T86K5Qa01nt1WjS\ngBycYNTL71aFk+DSmE7ZWLmHd1dUcMH4kpbBJbE0d+12lm2sYumGLmdRMsZ40F4fzXvAeyLypKp+\nFcMymR7O9w188sj4Ls/w9bYaxrgTOLsTVUUVkpLaBmnLrGMSkZc+mhp3PZoZIvKO7xb1kpkea9WW\nXfTLTWf8kNiknQn0zcMGA9DY3D0/lQ/71dsc95tZQZ9TGxBqEpCXQPMM8AUwFLgdWIuzlowxHbal\nupbPN1Rx0YRBcWk2A7jmxOGAk9CzO9pSXce67cH7mKxGYxKRl0BTqKqPAQ2q+p6qfg84PMrlMj3U\ns3PWoaqcN64kbmXIy3TS0HTXQOMvMAt1N62kmR7OS6Dx/TduFJEzRWQsztLJxnTY7FVbGT0wNy6j\nzXwyUpNJS0miqpsHmsamZobe1Go5J2Z8tjFOpTEmNC+ZDH/pJtT8Cc78mVzguqiWyvRI67bXMHfN\ndq47aUS8i0J+Zmq3rNFs81uG+urnPm3z/GfrbSSdSTxekmr+x71bCURlWWezb3hj6SYApsQwiWYo\ned000PzspSUt92d8timOJTHGu/YmbP5MVX8tIn8kSG4zVb06qiUzPc6/F63noOI8SuPYbOaTl5nK\nzpruF2h21NSH3eeuV5fxw2P3o2+vjBiUyJjw2uujWe7+nA8sCHLrMhE5TURWiEiZiNwY5Pl0EXne\nfX6OiJT6PXeTu32FiJwacFyyiHwqIv8JPKeJjw079/D5+irOOST+tRlwA003rNE0eejtf/SDNdzz\n2oqw+xkTK+1N2HzFvbtEVds2BneRiCQDDwInA+XAPBGZrqrL/Hb7PrBDVfcXkanAPcDFIjIamAqM\nAQYCb4nICFVtco+7BidQ5ka63KZzZq3YAsDRw2OfRDOYwpw0Pu8mmQH21DdR19jEm8s2s6vOW462\n1OT4DB03JhgvgwHuF5EBwIvAc6q6NELXngSUqepqABF5DpgC+AeaKcBt7v2XgAfEmXwxxS1LHbBG\nRMrc830sIiXAmcBdQKu1dEz8vDi/nJH9ejEqQfJw9clJZ9uuepqbNegM+0Qy6VdvdTjbdK47hNuY\nRBB2eLOqHg9MBiqAR0TkMxH5eQSuXQys83tc7m4Luo+qNuIMSCgMc+zvgJ8BzREoo4mAlZurWbRu\nJxdOiE9us2D65KTT2KwJPSCgqraB/3lpSaeWNEiQl9kYwNs8GlR1k6r+AbgCWATcEoFrB/tXCGyA\nDrVP0O0ichawRVXD9iGJyOUiMl9E5ldUVIQvrem0+978kpQk4RtjA79HxM/A/EwA1u2oiXNJQvvL\n+6t5fv66dvcJVRmrqWsK/oQxcRA20IjIASJym4gsBR4APiIyEzbLgUF+j0uADaH2EZEUIA/Y3s6x\nRwHniMha4DngBBH5W7CLq+ojqjpBVScUFcU3uWNP1tjUzDsrtnDeuGIKc9LjXZwW+/fNAaBsy644\nlyS03fXhg8WF4/f+G/z50nEt95/+5Cs2VvacpRBM9+alRvMEsAM4WVWPU9U/qeqWCFx7HjBcRIaK\nSBpO5/70gH2mA9Pc+xcA76iTc2M6MNUdlTYUGA7MVdWbVLVEVUvd872jqpdGoKymk5ZtrKK+sZkj\n4rDAWXv65zlDf7f6TYBMNLUN4QNNcW+nZvaNscUUBQxnPuL/nNy3u+oauevVZZ7OZ0w0eOmjORx4\nBGctmohx+1yuAt7AGSH2gqouFZE7ROQcd7fHgEK3s/964Eb32KXACzgDB14HrvQbcWYSyD8WlJOW\nksQJI/vFuyitZKY6a/fVeKg1xFpzs1JRXUdtQ/BuxuNG7K2BF7tNgDtr6slOD74e4R/fXsmjH6zh\nxTDNcMZEi5ems7Nx+mVedx8fKiKBNY9OUdUZqjpCVYep6l3utltUdbp7v1ZVL1TV/VV1km+Emvvc\nXe5xI1X1tSDnfldVz4pEOU3nVNU28I+F6zn9wP7kZSXWKKjkJCEjNSkhA83D769m4l1vUVYRvFnv\nz5eOb7k/ptgZwT+8Xy+y09oOIm3yG/Dwi5cjNWDUmI7xMrz5Npyhw+8CqOoi/4mTxoTy/Nx17Kpr\n5AdH7xfvogSVnZbCbo/zUmLpvS+dlulVIfqPMtOS6ZebzuaqOkb1z+U//300I/r1orq27Qi6ddtr\neO9LG+xi4stLoGlU1cpEGZZquoemZuXJj9Zy2NACDipJzFUss9KTE7JGk5rsNDT4T87cryib1RW7\nWx5Pv+poVrk1ngOLnde3IDuNn502km276nnswzUATH3kEzZV1caq6MYE5WUwwOci8k0gWUSGu7nP\nPopyuUw3N2/tdtbv3MO3Dh8S76KElJWaQk194tVoUoKMWa5raObp709i1g2TAeiXm8GRw1pnWRAR\n/mvy/vTP3TsoIDDIfPeJuZEvsDFheAk0/42T6qUOeBaoAq6NZqFM9/fqko1kpCZx4qi+8S5KSL0y\nUti+O3ySylhLSd77b3nqGGcQRW1DE8cML/K0js+4dpbIfneFNaOZ2PMy6qxGVW9W1YnuvJObVdXq\n4iak2oYmXv1sIyeO6kd2upfW2fgYOzifxesqqW9MrCQSM5dtbrl/1fHOstM/Os57P9f4Ib354Ge2\noodJHO0tE/AKQZYH8FHVc0I9Z/Ztby/fwvbd9Vw4IbEXYi3tk019UzPbd9e3zKuJt8CgV9QrnbV3\nn9nh85S482uMSQTt1Wh+A9wHrAH2AI+6t13A59EvmumOmpuV619YBMCE0oI4l6Z9BVlpAGzbnTiT\nNgNHjuV3cli4iHD/RYcEfc7m05hYa2+ZgPcAROROVT3W76lXROT9qJfMdEtvLd9MXWMzPzh6KDkJ\n3GwGzigtIKH6aar8EmheNKGEjNTgkzC9SEsJ/j3ypy8t4dBB+QzvlxiZtE3P52UwQJGItDQQuylf\nLDmYaUNVeeT91RTnZ3Lj6aPiXZywCnMSMND4ZZO++7yDu3SutOTQ/95nP/Bhl85tTEd4+cp5HfCu\niPhm5ZcCl0etRKbb+njVNuZ/tYPbzh7dauRUoirIdpJ8btuVOIHm7S/2phHs6jo5oWo0QMj0NsZE\nQ9hAo6qvi8hwwPcV9Qt3wTFjWpm1YgtpyUlcPHFwvIviSX5mKkkCO2riH2h21TXy6dc7KI/gsgXt\nBRqA8h01lPTOitj1jAnF63o0daq62L1ZkDFtbK6q5eVFGzhsvwIy0zrfrxBLSUlC76w0tiVA09l1\nzy/i24/NZZXf7P+uCtZ05h98Xv98U8SuZUx7Er99w3QLt7+ylG2767n2pBHxLkqHFGSnsT0Bms5W\nbq4GYGt15L7HBWt68x8+/ctXlyf0CqOm57BAY7qsorqON5du5ntHlTK+nVnpiaggOy2hBgPsimCS\nT1+YOaQkj/PGFvPYtAlt9lmzNXI1KGNC8RRoRKRYRI4UkWN9t2gXzHQfv3/7SxqblUsmdY++GX99\nctLZmkDzaCJZw/AlwlXg/osP5cQD+nHvBa1Hsp374Gx+8NQ8/rMkcHFbYyIn7GAAEbkHuBhnkTFf\nqlsFbC6NYeXmav72yddMGNKb/Ypy4l2cDku0Gk20XThhED99aUmrbW8t38Jby7dw1sED41Qq09N5\nGd58LjDSBgGYYJ6Z8zUAv7340DiXpHMKstPYWdNAQ1NzS3r+eLtk0mDOOKh/l89zwIBeHLV/ITee\ndkAESmVM53n5z1oNJNbyiCYhbKmu5R8Lyzn7kIEMKuiew2T7uJM24z3E2X+9p5+dOpJjhnd9TnR6\nSjLP/ODwkOsB/fzM1gGoLMRCa8Z0lZdAUwMsEpGHReQPvlu0C2YS303/+IzGJuXqE/aPd1E6zTdp\nM57NZ7vrGmls3jsarLP5zToqPSC9zUn3vxeT65p9j5dAMx24E2exswV+ty4TkdNEZIWIlInIjUGe\nTxeR593n5/gvIS0iN7nbV4jIqe62QSIyS0SWi8hSEbkmEuU0ba3Zupu3v9jCjycP69Y5s3xpaOKZ\nHWDMrW+wbvuelsfRXs32/osOYdoRQ8gIM6HTmEjxkhngKRFJA3wTJFaoapeHxohIMvAgcDJQDswT\nkemqusxvt+8DO1R1fxGZCtwDXCwio4GpOAuyDQTeEpERQCPwE1VdKCK9gAUiMjPgnKaLKmsauPQv\nc0hOEi6eOCjexemSwmxfBuf4BBrVkCtxRM1540o4b1yJjTQzMRP2K42ITAZW4gSFh4AvIzS8eRJQ\npqqrVbUeeA6YErDPFOAp9/5LwInifN2bAjznZixYA5QBk1R1o6ouBFDVamA5UByBsho/d766jPU7\n93Dr2aPpl5sY67h0VmGO23S2Kz5jXRqaYh9ofIItGW1MNHipO98HnKKqx7nLBZwK/DYC1y4G/BfG\nKKdtUGjZR1UbgUqg0MuxbjPbWGBOBMpqXEvKd/LSgnL275vDd44ojXdxusyX7yxeNZraxqbwO0VJ\nU5C8mjt21/OLf39ObUP8ymV6Hi+BJlVVV/geqOqXRGYUWrCvU4Ff70Lt0+6xIpID/AO4VlWrgl5c\n5HIRmS8i8ysqbB11L1SVX766nLSUJJ747sR4Fyci4p3vLJ4f6P4DEHzuef0Lnv7kK15ZbM1qJnK8\nBJr5IvKYiEx2b48SmcEA5YB/A38JEPjubtlHRFKAPGB7e8eKSCpOkHlGVf8Z6uKq+oiqTlDVCUVF\ntryOFy/OL2fumu1cOXn/bjucOZjCnDS2xanprC4gXf83xsaupTc7rW0X7YbKWiD6AxLMvsVLoPkx\nsBS4GrgGJ0PAFRG49jxguIgMdQcbTMUZ4eZvOjDNvX8B8I46vafTganuqLShwHBgrtt/8xiwXFXv\nj0AZjWvO6m3c/O/PmFRawFXdeDhzMPHMDuBfo8lKS47pxNcTD+jLPecfxJDCvV8a3v/Sqd1b/42J\nJC+jzuqA+91bxKhqo4hcBbwBJAOPq+pSEbkDmK+q03GCxtMiUoZTk5nqHrtURF7ACXqNwJWq2iQi\nRwPfBj4TkUXupf5XVWdEsuz7mt11jVz17Kf0y83gT5eOI7mHfQgV5qSzfGPQFtao2+MXaHIzYjsv\nWkS4eOJg5q7ZwVfbWq+D09VF14zxFzLQiMgLqnqRiHxG274TVLVr68w655gBzAjYdovf/VrgwhDH\n3gXcFbDtQ4L335gueGBWGRXVdfzjx0e2jNLqSQqz0+I2j8Z/pcvmOAx1Brjz3DHsrKlvtbpnsjWd\nmQhqr0bjm+x4ViwKYhJPc7Nyzxtf8PB7q7loQkm3WwLAq4LsNCr3xCffmf/6MPEKNFlpKTz23Yn8\nc2E517+wOC5lMD1byP8qVd3o3v0vVf3K/wb8V2yKZ+Lp73O/5uH3VnPSAX2589wD412cqPHV0mKd\n72xLVS3fe3JeTK/Znpz0vd87r/z7QnbXNdLcHL95Pqbn8PL17eQg206PdEFMYqnc08B9b67gsKEF\nPPqdCaSndI/lmTujIMtNrLk7tqtN/uTFxdS7k1nGDc7nycsmxfT6gZoCgsqYW99gv/91WrZ/+uJi\nfvPGimCHGRNWe300P8apuewnIv4LWPQCZke7YCZ+1u/cw0V//pjKPQ384qzRPX6oa1aaE0T3xHhO\ni/+girvPP5gRcc4Zd9LofiGfe3FBOQCXHj6E/nndOxuEib32ajR/B87GGUp8tt9tvKpeGoOymTio\nrGlg2uNzWb9zD/decAgHFgdPMd+TZLhZjPfUxzbQ+AeWRFgLJzU5iU9/0bYB4+mP17bcf2H+ujbP\nGxNOyBqNqlbipHy5BEBE+gIZQI6I5Kjq17EpoomV2oYmzvzjB5Tv2MOvzz+Y88eXxLtIMZHp1mhi\nPUvffyBAanJi1BpTg2R0vuf1vU1mgwoyY1kc00N4Sap5toisBNYA7wFrgdeiXC4TY1uqazn1d+9T\nvmMP95x/EBd186zMHZGZGp+ms7pWgSb+NRqAtCDl2FXX2HI/qYc3o5ro8PLu/iVwOPClqg4FTsT6\naHqUFZuqOeaeWXy1rYbrTx7BxRMHx7tIMeULNDUxbjqr80uomSiBJlzNyprOTGd4eXc3qOo2IElE\nklR1FtA9F4g3bXy4civnPTSbusZmrj95BFefODzeRYq5jDTn3yDWNZpEbDoLN/Bjdtm2GJXE9CRh\nU9AAO91syO8Dz4jIFpy0L6Yba2xq5onZa7n3jRUM7ZPNQ5eOY1hRTryLFRdZbnLJPfWxfVvXJ2DT\nmTHR4OXdPQXYA1wHvA6swhl9ZrqplZurmXjXW9w1YzlHDCvk+R8dvs8GGYCs1GRy0lNaLaccbZU1\nDby5bHPL4+4UaG540bIHmI4J++5W1d2q2gRkAa8AfyNI7jPTPTw4q4zzHvqIHTUN3Hj6KJ68bCL5\n7oTFfVVSkjBmYC6fb6iM2TV/9/aXrR4neqLS40bsXUrjJXdOjTFehW06E5EfAXfg1GqacZJWKrBf\ndItmIkVVmV22jV/NWM6yjVUcN6KIm888IO4TBBPJoIIsPly5NWbXq9wT2ywEXfXAN8dy0G1vtjxu\nbGompRvVwkx8eemjuQEYo6qx+y80EfPWss384K/zAeiVnsKPJw/jhlNGJvw36Fjrn5tBxa66qH+A\nfr2thmPvnUWfnMStRZb0zqR8x95mxILsNHoFLGFw+yvLenT+OxNZXgLNKqAm7F4mYagqLy0o5+f/\n/rxlrsa5hw7k9ikHkpcZ2zVPuot+uek0NSvbd9fTNzd6KVbeWu70y2zdVc9xI4p44Jtj47boWijv\n/GQyivLNR+cwsbSAG08f1Waf1z7faIHGeOYl0NwEfCQic4CW9W5V9eqolcp0yrZddTz50Vr+9en6\nlm+kJx3Ql2tPGrFPpJLpilw1utboAAAgAElEQVQ3AFfVNkQ10PiPHu6dlUqvjNQ2tYV4S3OzA/zj\nx0eG3GfrrnqmL97AgQNz2S9gIElDUzN1jc2tskGbfZuXd8LDwDvAZzh9NCaB1DY08Z8lG3n647Us\nLt/bmX3GQf35xVmjGZBnKUO88NX0ot134t9g2d0GYRw5rJCPVu2dR3P1s58C8JOTR/DcvHXMvvEE\nAK58ZiFvLtvM2rvPjEs5TeLxEmgaVfX6qJfEeKKqrKrYxcxlW1i6oZK3lm9uWaVxcEEWv/rGQYwd\nnE+2fZvskFgFGv8lkguzu1egufH0UZzzQNukIPfNdEbQ+fq3/IdtGwPeAs0sEbkcZ2izf9PZ9qiV\nygDOP+6mqlrW79jD7FXbmLdmO6u37mJzlfNnKMxO47xxJUws7c1JB/RLuCaY7iRWgcZ/Ec1+UWyi\niwbfcgqh/P7tlUweuXcYtKq2yjSwfuceBBiYb7XsfY2XQPNN9+dNftsiMrxZRE4Dfg8kA39R1bsD\nnk8H/gqMB7YBF6vqWve5m4DvA03A1ar6hpdzJqL6xma2VNeyuaqW5Rur+WJTFZsq6/ho1dZW+bcK\ns9MYN6Q3R+/fh6OH92FoYXarb8im84p6OatsbqmqC7Nn1/inuembmx7Va0VaZlr7Hxd/fKeMP75T\n1vK4Ylcdk+56mz99axynHzSAo+5+B8Ca1PZBYQONm0gz4kQkGXgQZwXPcmCeiExX1WV+u30f2KGq\n+4vIVOAe4GIRGQ1MBcYAA4G3RGSEe0y4c8aEbwRTdW0DFdV1VNc2snVXHVt31bGpqpZNlbXuT2db\noF4ZKRw3ooje2Wkcs38fDtuvkIJu1tTSnfTKSCU3I4X1O6ObHcD3xSE9JYlR/XOjeq1Iy0rt2Cqr\nKzZVA/Dw+6s5/aAB0SiS6SbaW2HzBFV9R0TOC/a8qv6zi9eeBJSp6mr3es/hpLvxDwpTgNvc+y8B\nD4hTF58CPKeqdcAaESlzz4eHc0bME7PXMGtFBQ2NzWzdVUdNfRO1DU3scW8aIn9CflYq/XMz6J+X\nwYED8+ifl0H/3Az65WXQJzudgpw0iq15IeYG5meyYWdt1M6/fXc9f3h7JQCzbzyBPjndrUbTsUDT\n6C4NraH+Ecw+o70azXE4o82C5TVToKuBphjwzzleDhwWah9VbRSRSqDQ3f5JwLHF7v1w54yYxial\nck8DyQLDinLISk8mMzWZrLRkMtNSKMhKJTs9hQF5meRkpFDUK53C7LSWFR1NYumdlUblnujNaZld\ntnfOc3cLMuDUwjri8Q/XANAcEGd+NWM51540vCWZqT9VZdoT87jsyFKOH9W302U1iaW9FTZvdX9e\nFqVrB+tcCPzqE2qfUNuD/ScE/TrlDnC4HGDw4M6tv/LDY/fjh8daJp6eIi8zldVbd0Xt/Au+2hG1\nc8eCiDCptIC5a72NA/rATenTHFCjeeT91WSlJXPtSSPaHFPX2Mz7X1bw8aqtrLzrjK4XOkFt313P\nuDtn8tC3xnHGPtCs6GWFzadFJM/v8RAReTsC1y4H/JdxLAE2hNpHRFKAPGB7O8d6OScAqvqIqk5Q\n1QlFRUXBdjH7mPysVHbWRG/U2ZLynVE7d6xoJ/LpLt1Q1WqRN4DmwGqOy7ecdrh1cbyorGlg//+d\nwQcrK1q2rdm6m9c+29ip86lqxPrwyrY4X2h8tb6ezktd+ENgjoicISI/BGYCv4vAtecBw0VkqIik\n4XTuTw/YZzowzb1/AfCOOg2+04GpIpIuIkOB4cBcj+c0Jqi8rNSoDm/eXed8iH5jbHGYPRNXZ7tb\n/m/GF60e+0ZLrtte0yqZqW9OWHuDKW+bvpT73lwR9pqfra+ksVl5aNaqlm2n/PY9fvzMwo4UvcXr\nn2/iqLvf4b0vK8LvHIbv92vaR/qvvIw6e1hElgKzgK3AWFXd1NULu30uVwFv4AxFflxVl4rIHcB8\nVZ0OPAY87Xb2b8cJHLj7vYDTyd8IXOkuZUCwc3a1rGbfkJeZSl1jM7UNTVHpR9tV18j540q476JD\nIn7uWPGNfExOEprcWsmD3xzHlX9v/8P7k9WtV+ZMcT9pJ//mXZqalecvP5z7Z37ZMuAgqZ0azZMf\nrQXgJ6eMbPeaDc1O0Er161tqaNo7QCFcrWn9zj0MzMvgw7KtJCcJ29ycdC8vWt9q2QQvauob+WT1\nNk4Y1Q/YuyxEqJpdJO2qa6SmvpG+veI3b8vLMgHfBn4BfAc4GJghIpepapdXP1LVGcCMgG23+N2v\nBS4McexdwF1ezmmMF/mZzofozpoG+udFPtBU1zbQK6N7Z2z49QUHc8ziDZw3roSK6jr65WawqSr8\nSL0v3KHOPslJzoe/L1hd/MgnrZ4PDAGPvr+avrnpFGbvHURxxyvLmDyyiGNDfOg3ukElNUkovfFV\nfnrq3sDU2Kytls+eXbaVJBGOGFYIwMKvd3DeQx/x6/MP5mf/WALA7y4+tNV5axuaSE1O8pQJ/ZaX\nl/LSgnJev/YYRvXPbWmA7Eic+fTrHQwqyOrwQJIz//ABX22riev8JS/v+vOBo1V1C/CsiPwLeBIY\nG82CGRNr/tkB+udF9tufqrKrrrHbJ5rMz0rj20eUArSkOSrsxJIHryzewMxloRtGdtc3sae+qaWG\nc9eM5W32eXz2Gh6fvSbkB2hDU+vUjPe+sbe5rb6xudWqpt/6yxwA/vq9SRw7oqilD8W/mcw3XFvE\n+XuO+sXrXDi+hHsvDF5DraxpYOydb/LEZZNaFovbWFlLU7Ny3kMfAXsDrRffcI8JFzCWb6yiMDuN\nvrkZbN1Vx1fb4p9838sKm+e6Qcb3eC5RHDJsTLzkZzmBZmdN5Ic472loolkhp5vXaILpFRA8Swuz\nuGPKmHaPWbaxioVftz844kd/W8C3H5vDjk4uo+ALNPVNbXMB+5bPCPSdx+dSeuOrfLVtN0Cr2lqj\ne56qPQ2c8tv3AXhxQTmL1u1sec5fWUU1zQq/8Qtwlz0xj9/OXNny2BdoGpqaPQ80CDcv6fTff8Ax\nv54FwEV//tjTOaPNy6izEhH5l4hUiMhmEfkHYAPcTY8TzXxnvtQ23b1GE0xgX8eQwmy+49Z6uuL9\nLyv4YOVWZoZJ0rml2gkGlTUNNDQ189KCcmobmlr6Y+rdoOLfxFXX2MRbyzYzZ/W2oB/wS9xM6L5z\nA9z4z88AmLWigpVb9g6DP/fB2S2JRf2lJTu1scD1hvwDhW/o989eWsJRd7/TMuoukH/tbOuu4IF3\n++76loDnC6Srt+4Oum+seXnXPwH8nb19JZe6206OVqGMiQdfoNkZ4UBT29DE5N+8C9Dt+2hC+fKX\npzPi568B8IepkW1V9/WRhDLprrd587pjOeW379Pf7TO64cXFDCnMAvY2efkPYKhvbG5ZeTYY3xyg\nddu91TK+2FgFOHOlxgzMZd32mpYh3dt2t04x5T+vqEmVE+57l9UVTkD4z5KNXDC+hN+99SXTF23g\nnRsmA7QKQFuqa1ty84ETNNdt38NJ97/H5e3M63tpQTkXjC/x9PtEmpd3fZGqPuH3+EkRuTZaBTIm\nXvLcprOqCAeaT/2aiHpijQb2LpYGe1/HYC6ZNJhn534d8ev7+lj8m7p8fRP+wcUnVNNZZ320ahu3\nTV/aMiLOn2/Ito9/w9eGnXtaPX/Di4v5+5yvWpoVm5qVHzw1r1Ut5sOVWxkzMI/KmgZeWljOBysr\neHeF05f05Oy91//JC63Ha93w4uKEDjRbReRS4Fn38SU4mZSN6VF6paeQnCQRn7RZU9/Ycr+nBhqv\nrjt5eFQCTUV16Kzbi9a17Quqa4hsoKlrbA4aZILx7/8PDEJAq76rlxasY9aK1vN2/u+1LzjtwP4c\nd++7bY7174/6x8JyT+WJBS8TNr8HXARsAjbiTJz8XjQLZUw8iAh5mansiPBgAP9v0j1xMIDPgp+f\nxEfuKpuh+A9PDuU3IUZxRVJ9U/C+kFh4vwMTPv2XXfD31vItQbeHM/aON3lzaZenQXZYu4HGTeV/\nvqqeo6pFqtrXHYX2VYzKZ0xM9e2VzmYP80I6wjdZMSc9hUEFWRE9dyIpzElvtaiZb6G07xwxpGWb\nlzknscic8M+F66N+jY6YcujAoNuDLSEC8OqSoJm1wtpR08DlTy9g+cYqFny1vc0Q8Ghp9+uVqjaJ\nyBTgtzEpjTFxVpyfSfmOyK1JM2f1Np762Ple9vq1x5C7D62CmpWWQk19E1ccN4z8zFQme8jGPHXi\nIE/BqKuemeO9+a4gO63NyLFIm3ZkKS8vahs8gjWtAWGHhodz+u8/AOCHxwzl5jNHd+lcXnhpOpst\nIg+IyDEiMs53i3rJjImD4t6ZbIjg4mf+Q2fTU/at5SGy053ft66xmetPGcm4wb3b3X9EvxzuPv/g\nWBStjW8fPoSLJgTvKA+15PbE0vZ/n47Yv29OxM4VaPSAXE4Z3S/oc49+sIY99dFvRvQSaI7EWcny\nDuA+9/abaBbKmHgpzs+kqraR6trIDAjw7/hN6+B6Lt3dQ98axzmHDGRQ79aL+PVKT+G8ccV8+D/H\nt9re3iR5/+a3aPj2EUMo7ZPdZvuhg/LZVRf8vRDJv2d2WgrfPbK0w8fdctZo/n3lUe3us3xTFQ9/\ne3zI599f2fUkoeF4yQxwfJBb+z1+xnRTxe6HYqTSwfvPmejowmHd3ZiBefzhkrGkJLf+vT+7/VTu\nv+jQNkke25vxnp/prclx2hFDePO6Y4M+d9bBe9d9ueyoUkb179XyeP+inJbM1D88Zu/q9Y9/dyLN\nIboxUpO9/T0Hh+iXy/cbBp6cJK3mxoRy3UkjuM5vHZ9pR5Zy6KB8yu46vdV+mX5JYVXbX3bh2OHR\nXybFS2aAQhH5g4gsFJEFIvJ7ESmMesmMiYOCLCdv147dkanR+H94pnn8YNpXpAT0xQSLMyP7OcGg\nIUxOsAOLc7l4wiBun3IgI/r1CrqPLxnmr75xELeePYbXr90bkJKSpOVv5R8YU5KFxhCRJjvICqHB\nvP+z44NuP29s66Y6L3nPLppYQv+8vQHJ158VGMwH+OXq88WY/z1jVNBzdnSJ7s7w8s5/DqjASa55\ngXv/+WgWyph48SWK9J/70hX+H55JMejk7k4CX4/AlThh7+TPsYPyOdNvJco7poxpVSN55aqjufv8\ng9q9nm+m/oh+wftDJpYWAHDksL3fo1OTkgg1MOuOKWP47xP2b3n84hVHMGZgbtB9/csK8NT3JpHi\nZo8+dYzTf+JbabNfbuiaTWZqcst7KlSfEsAfLtmbneH9nzqB7vJjhzH/5yeFPCaavASaAlW9U1XX\nuLdfAvnRLpgx8eALNLvqIhNoGmOw3khPEeyV+vbhTt/MQSV5PPitvWOQvnNEaasaiYi02zw0fkhv\nfn7WaK4/eQTjh+ztxL/t7NH8/YdOjuDD9itk6e2ncoxfU1JKsnD+uLbDrd/+yXEU5qS3WhNnYmkB\nr159DGceNKBNYJlx9TE8+8PDWx4fN6KopTZycInzcbp/3xzW3n0m+/UJPTAgIzWZ40f1JTcjhe8e\nObTVc/59XoN6Z7HktlP42/cPazWkPl7Nt16uOktEpopIknu7CHg12gUzJh58M/d9q2F2Vagkiaat\nYE1nZx8ykLV3n8mAvMy2T7bjsWkTACd/3eJbT+HvPzyMYUU5XH3i8FYB6btHDeXIYX1aHvu+aExz\nBx+kJAn/c1rbJqdhRa2DwX5+Awke/Na4VkEQnNpbVkATVaobaALnsgSr2fmkpyTRLzeDJbedyuiA\n2lNJb7+AkppEbkYqRw/vE3B8fEY+emlk/BFwPfA393ESsFtErgdUVYPXFY3phrLcIbm7I1SjiebS\n0D1Nex+wHTV5ZF8OGJDL9SePaEmW2hG3nj2Gm88c7daU2t93+R2nkeThK3tgQtXAxd982nsdwq0K\nOm5wPgu/3hmyP9B/sbdY8rKUc/CeNWN6IF8HbySazmobmnhh/roun2df4f/5+tNTR5IbJF3P+z89\nvqWvpT3JScJr1xzT6bIkJQlpHvvUvHamD3VrPSe4E1d9fTSBzavBWlvzMlN55gfhlwF74UdHsLm6\nLmR/oC9QTRpaQEl+Jlf59TFFk6dhEyJyHnA0TjPqB6r676iWypg4SU4S8rNS26R274xlG6vY7K5D\nc2CxVfzD8R+hd+XxwT8ABxe2Hir82W2nRLVMXfHZbae06ncSEZbcdkpLP4lv1F1gjSbY6LNLDx/M\ngcV5Ya+ZkpxEcX77zYxf3Hma5yWoIyVsoBGRh4D92Zu9+QoROVlVr+zsRUWkAGfkWimwFrhIVXcE\n2W8a8HP34S9V9Sl3+3ic5aQzgRnANaqqInIvcDZQD6wCLlPVruVqMPucgXmZbNjZ9XxnKzdXA/Da\nNcewX1HbyYCmtWM6MZ+jV4xT+gwqyGS9xxRFwcrmn4IoOUQfjS/g3n7OGAqy0zh2RFGbVUy7IiM1\n9v00Xkp/HHCgur+9iDwFfNbF694IvK2qd4vIje7j//HfwQ1GtwITcGpSC0RkuhuQ/gRcDnyCE2hO\nA14DZgI3qWqjiNwD3BR4XmPCGZifEZF8Z5sqndrMsKKcfS4rQEccOayQ/zvvoA53+MfDrJ9MDjo6\nrjN8NZrALhlfheaQQfkcOqhnDPD1EmhWAIMBX8bmQUD7S96FNwWY7N5/CniXtgHhVGCmqm4HEJGZ\nwGki8i6Qq6ofu9v/CpwLvKaqb/od/wnOvB9jOqRPTjqL3aV8u6JiVy29s1ItyLRj5V2nkyQS02ac\nznj3hslkpSe3mRjZFRdPHMyyjVVcc+LwVtt9TWeBE1q7My+BphBYLiJz3ccTgU9EZDqAqp7Tiev2\nU9WN7vEbRSRYWtdiwL8ntdzdVuzeD9we6HvYxFLTCQXZaezYXY+qhh3l056K6jpPaUX2ZV7TuMRb\nsDxoXZWZlsyvL2i79o5v1FlSF957icZLoLmlMycWkbeA/kGeutnrKYJs03a2+1/7ZqAReKad8l2O\n0/zG4MGDPRbJ7AsKstNobFbKtuxieIh0Jl5YoDGd4Qs0iV7L6wgvw5vf838sIkcB3ww3GEBVQ+Y6\nEJHNIjLArc0MAIItF1fO3uY1gBKcJrZy977/9paFHNwBBGcBJ6qGHpCuqo8AjwBMmDDBpm+bFr6J\ncH/75Ctun3Jgp89TsauO8WFS4xsTyNdH04PijKfMAIjIoSLyaxFZC/wSWN7F604Hprn3pwEvB9nn\nDeAUEektIr2BU4A33Ca3ahE5XJx2je/4jheR03D6es5R1ZoultHso44c1ocJQ3qzfFN1p8+hqlaj\nMZ3yx0vGcu6hA1vm3fQEIWs0IjICmApcAmzD6e8QVQ2eirRj7gZeEJHvA18DF7rXnABcoao/UNXt\nInInMM895g7fwADgx+wd3vyaewN4AEgHZrpt65+o6hURKK/Zx+RlprKxsvNDnK/42wJqG5ot0JgO\nO2BALr+bOjb8jt1Ie01nXwAfAGerahmAiFwXiYuq6jbgxCDb5wM/8Hv8OPB4iP3atGmoamymuZoe\nLycjhd0Vnc8O8MbSzQAUZFugMaa9prPzgU04STUfFZETCd4Rb0yPk52ewq7ayOQ7M2ZfFzLQqOq/\nVPViYBROJ/x1QD8R+ZOIJG7eB2MiICc9pdP5zvzHoBw2tCBSRTKm2/KylPNuVX1GVc/CGeG1CGcm\nvzE9VnZaCnWNzW3Sg3ixeutuAH5+5gGt1gIxZl/VodlSqrpdVR9W1ROiVSBjEsFQNzfZkvKOp8r7\n1qNzgO4zGdGYaLP/BGOCOKTEyZS7ZmvnR8kfVBI+264x+4LIpQQ1pgfJctelqanveD/N4MIsBuRn\nMM4maxoDWI3GmKCy3ZU2a+o7thRzfWMzX26ubnfdd2P2NRZojAkiI6Vzgebr7bvZWdPA0cMLo1Es\nY7olCzTGBJGUJGSmJlPTgSHOX23bzbTHnUQW+Vlp0SqaMd2OBRpjQtjT0MRfPlzjef87/7OM9Tud\nBdOCrXdvzL7KAo0xEZLut0RurJcYNiaRWaAxJoRxg51ldOsbvU3azPQLNDkRXOPdmO7OAo0xIZx1\n8EDA+xBn/0DT2/pojGlhgcaYEHxDnMt37PG0v2/l3VW/OoPMtOT2dzZmH2KBxpgQMt1Jm1MenO1p\n/111jRTnZ/aoJXiNiQQLNMaEsNsd2tzU7G2l76o9jfSy0WbGtGGBxpgQ+udmdGj/iupaW1HTmCAs\n0BgTwvGj+pKZmszR+/cJu+9n5ZUsLq+kXweDkzH7Ags0xrRjQmlvTwugnf3AhwAMzM+MdpGM6Xbi\nEmhEpEBEZorISvdn0DS3IjLN3WeliEzz2z5eRD4TkTIR+YOISMBxN4iIikj4r6LGtCMnPYVF63ZS\nXdsQcp/ahr350A7fz1bUNCZQvGo0NwJvq+pw4G2CrNgpIgXArcBhwCTgVr+A9CfgcmC4ezvN77hB\nwMnA19H8Bcy+oTDHmQ9z35tfhtznJy8ubrl/+FBLpmlMoHgFminAU+79p4Bzg+xzKjDTXdVzBzAT\nOE1EBgC5qvqxOouz/zXg+N8CPwO8DRUyph3/fcJwAOoaQ2dxfnPpJgB+cPRQkmxoszFtxCvQ9FPV\njQDuz75B9ikG1vk9Lne3Fbv3A7cjIucA61V1McZEQL/cDPrlprc7xLmhyXmutp1gZMy+LGqBRkTe\nEpHPg9ymeD1FkG0aaruIZAE3A7d4LN/lIjJfROZXVFR4LJLZF+2ua+KF+eU4FejWVlXsarmfk26J\nNI0JJmqzy1T1pFDPichmERmgqhvdprAtQXYrByb7PS4B3nW3lwRs3wAMA4YCi92xASXAQhGZpKqb\ngpTvEeARgAkTJlgzmwnJN+psU1UtA/JajyrbXFnbcv+aE4fHtFzGdBfxajqbDvhGkU0DXg6yzxvA\nKSLS2x0EcArwhtvUVi0ih7ujzb4DvKyqn6lqX1UtVdVSnIA0LliQMaYzgmVx3u23AqflNzMmuHgF\nmruBk0VkJc4IsbsBRGSCiPwFQFW3A3cC89zbHe42gB8DfwHKgFXAa7EtvtmXPPDNsQDUNrQNNFV7\nQg97NsY44pKYSVW3AScG2T4f+IHf48eBx0Psd2CYa5R2uaDGABkpTk0l2MizjZXeMjsbsy+zzADG\nhJHhrjMTrEbz9faaWBfHmG7HAo0xYaSnOv8mW6prWe03ygxgU1VdPIpkTLdiOc2NCcPXdHbV3z8F\n4M+Xjuf4UUU8/N5qtlTVtneoMQYLNMaElZHauuJ/xd8WcOH4El5cUB5yH2PMXvbfYUwYEmSK8Jbq\nvU1mp47px7ybQ04bM2afZ4HGmDBKemcxtE92q23vfbk3m0Rxfha9MiwrgDGhWKAxJoyM1GRm3TCZ\nK44bFvT5ZPsvMqZd9i9ijEclvYMvanbZUUNjXBJjuhcLNMZ4lJ7S9t/lmOF9bFVNY8KwUWfGeFRV\n23pJ599PPZQphxbHqTTGdB9WozHGI/8lmwGOG1EUp5IY071YjcYYj6YdWUrlngYumTSYT1ZvIz8r\nLd5FMqZbsEBjjEc56Sn87xkHALQZ7myMCc2azowxxkSVBRpjjDFRZYHGGGNMVFmgMcYYE1UWaIwx\nxkSVBRpjjDFRZYHGGGNMVFmgMcYYE1WiqvEuQ9yJSAXwVScP7wNsjWBxIsXK1TFWro6xcnVMopYL\nula2IaoaNheTBZouEpH5qjoh3uUIZOXqGCtXx1i5OiZRywWxKZs1nRljjIkqCzTGGGOiygJN1z0S\n7wKEYOXqGCtXx1i5OiZRywUxKJv10RhjjIkqq9EYY4yJKgs0HojIhSKyVESaRWRCwHM3iUiZiKwQ\nkVNDHD9UROaIyEoReV5EIr5ilnveRe5trYgsCrHfWhH5zN1vfqTLEeR6t4nIer+ynRFiv9Pc17BM\nRG6MQbnuFZEvRGSJiPxLRPJD7BeT1yvc7y8i6e7fuMx9L5VGqyx+1xwkIrNEZLn7/r8myD6TRaTS\n7+97S7TL5V633b+LOP7gvl5LRGRcDMo00u91WCQiVSJybcA+MXu9RORxEdkiIp/7bSsQkZnuZ9FM\nEekd4thp7j4rRWRalwujqnYLcwMOAEYC7wIT/LaPBhYD6cBQYBWQHOT4F4Cp7v0/Az+OcnnvA24J\n8dxaoE8MX7vbgBvC7JPsvnb7AWnuazo6yuU6BUhx798D3BOv18vL7w/8F/Bn9/5U4PkY/O0GAOPc\n+72AL4OUazLwn1i9n7z+XYAzgNcAAQ4H5sS4fMnAJpx5JnF5vYBjgXHA537bfg3c6N6/Mdj7HigA\nVrs/e7v3e3elLFaj8UBVl6vqiiBPTQGeU9U6VV0DlAGT/HcQEQFOAF5yNz0FnButsrrXuwh4NlrX\niIJJQJmqrlbVeuA5nNc2alT1TVVtdB9+ApRE83phePn9p+C8d8B5L53o/q2jRlU3qupC9341sBwo\njuY1I2gK8Fd1fALki8iAGF7/RGCVqnZ2IniXqer7wPaAzf7vo1CfRacCM1V1u6ruAGYCp3WlLBZo\nuqYYWOf3uJy2/4iFwE6/D7Vg+0TSMcBmVV0Z4nkF3hSRBSJyeRTL4e8qt/ni8RBVdS+vYzR9D+fb\nbzCxeL28/P4t+7jvpUqc91ZMuE11Y4E5QZ4+QkQWi8hrIjImRkUK93eJ93tqKqG/7MXj9fLpp6ob\nwfkiAfQNsk/EX7uUrhzck4jIW0D/IE/drKovhzosyLbAYXxe9vHEYxkvof3azFGqukFE+gIzReQL\n95tPp7VXLuBPwJ04v/OdOM163ws8RZBjuzwc0svrJSI3A43AMyFOE/HXK1hRg2yL2vuoo0QkB/gH\ncK2qVgU8vRCneWiX2//2b2B4DIoV7u8Sz9crDTgHuCnI0/F6vToi4q+dBRqXqp7UicPKgUF+j0uA\nDQH7bMWptqe430SD7RORMopICnAeML6dc2xwf24RkX/hNNt06YPT62snIo8C/wnylJfXMeLlcjs5\nzwJOVLdxOsg5Iv56BQwimAEAAAOZSURBVOHl9/ftU+7+nfNo2ywScSKSihNknlHVfwY+7x94VHWG\niDwkIn1UNap5vTz8XaLynvLodGChqm4OfCJer5efzSIyQFU3uk2JW4LsU47Tl+RTgtM/3WnWdNY1\n04Gp7oigoTjfTOb67+B+gM0CLnA3TQNC1ZC66iTgC1UtD/akiGSLSC/ffZwO8c+D7RspAe3i3whx\nvXnAcHFG56XhNDtMj3K5TgP+BzhHVWtC7BOr18vL7z8d570DznvpnVDBMVLcPqDHgOWqen+Iffr7\n+opEZBLOZ8q2KJfLy99lOvAdd/TZ4UClr8koBkK2KsTj9Qrg/z4K9Vn0BnCKiPR2m7pPcbd1XixG\nP3T3G84HZDlQB2wG3vB77macEUMrgNP9ts8ABrr398MJQGXAi0B6lMr5JHBFwLaBwAy/cix2b0tx\nmpCi/do9DXwGLHHf5AMCy+U+PgNnVNOqGJWrDKcdepF7+3NguWL5egX7/YE7cAIhQIb73ilz30v7\nxeA1OhqnyWSJ3+t0BnCF730GXOW+NotxBlUcGYNyBf27BJRLgAfd1/Mz/EaLRrlsWTiBI89vW1xe\nL5xgtxFocD+/vo/Tr/c2sNL9WeDuOwH4i9+x33Pfa2XAZV0ti2UGMMYYE1XWdGaMMSaqLNAYY4yJ\nKgs0xhhjosoCjTHGmKiyQGOMMSaqLNAYY4yJKgs0xhhjosoCjTFxIiK74l0GY2LBAo0xxpioskBj\nTJyJyPUi8rl7u9Zv+y/EWQV0pog8KyI3+D03RkTeEpEv3f3+KCIT4/MbGNM+y95sTByJyHjgMuAw\nnPxcc0TkPZwVGs/HWQMmBSe9/AL3GF/eswtxVj/8AligqvNi/gsY44EFGmPi62jgX6q6G0BE/omz\neF0S8LKq7nG3v+J3zEnAp6q61H0uDWedH2MSkjWdGRNfoZZjbm+Z5rE4NRxEZCCwS1VnR7pgxkSK\nBRpj4ut94FwRyXLXVvkG8AHwIXC2iGS4K1ye6XdMHc5iVAD/B6TFssDGdJQ1nRkTR6q6UESeZO+C\neX9R1U8BRGQ6zrolXwHzgUp3n78DL4vICuBhIF1Efqeq12JMArL1aIxJUCKSo87a8lk4NZ/LVXVh\nvMtlTEdZjcaYxPWIiIzGWV3zKQsypruyGo0xxpiossEAxhhjosoCjTHGmKiyQGOMMSaqLNAYY4yJ\nKgs0xhhjosoCjTHGmKiyQGOMMSaqLNAYY4yJqv8HIsbIhQRB2r0AAAAASUVORK5CYII=\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x28a0077b198>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.plot(log_alpha_grid, mdkl_sampled - mdkl_approx(log_alpha_grid, k1, k2, k3))\n",
    "_ = plt.axes().set_xlabel(r'$\\log\\alpha$') \n",
    "_ = plt.axes().set_ylabel(r'Approximation deviation') "
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.6.3"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}
